
R version 4.1.2 (2021-11-01) -- "Bird Hippie"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[Previously saved workspace restored]

> ## Code to replicate Table 2 of Abrajano, Elmendorf, and Quinn. "Measuring
> ## Perceived Skin Color: Spillover Effects and Likert-type Scales". JOP.
> 
> 
> library(clubSandwich)
Registered S3 method overwritten by 'clubSandwich':
  method    from    
  bread.mlm sandwich
> 
> 
> sink(file="AbrajanoElmendorfQuinn-Table2-output.txt")
> 
> 
> 
> ## JUST BLACK PHOTOS
> ## Part 1
> ## E[y | i %in% I^W] > E[y | i %in% I^B]
> ##
> ## where the expectation is over all respondents assigned to M1 or M2,
> ## all positions in S, and all photos in J^B
> 
> mydata <- read.csv("ScaleRaceSpring2017clean.csv")
> 
> ## just black respondents
> mydata.b <- mydata[mydata$R.race == "Black",]
> mydata.b <- mydata.b[!is.na(mydata.b$R.race),]
> ## just white respondents
> mydata.w <- mydata[mydata$R.race == "White",]
> mydata.w <- mydata.w[!is.na(mydata.w$R.race),]
> 
> ## keep all photos but nothing else
> inds <- 1:24
> keep.vars <- c(paste("X", inds, "_Q339", sep=""),
+                paste("X", inds, "_Q329", sep=""),
+                paste("X", inds, "a.id", sep=""))
> 
> mydata.b <- mydata.b[, keep.vars]
> mydata.w <- mydata.w[, keep.vars]
> ## convert factor to character
> for (i in 49:72){
+     mydata.b[,i] <- as.character(mydata.b[,i])
+     mydata.w[,i] <- as.character(mydata.w[,i])    
+ }
> ## get rid of white photos
> for (i in 1:24){
+     i2 <- i+24
+     i3 <- i+48
+     indic.b <- grepl("mw", mydata.b[,i3])
+     indic.w <- grepl("mw", mydata.w[,i3])
+     mydata.b[indic.b, i] <- NA
+     mydata.w[indic.w, i] <- NA    
+     mydata.b[indic.b, i2] <- NA
+     mydata.w[indic.w, i2] <- NA    
+ }
> mydata.b <- mydata.b[,1:48]
> mydata.w <- mydata.w[,1:48]
> 
> ## reshape into long format
> mydata.b.long <- reshape(mydata.b, direction="long", varying=1:48, v.names="MM")
> mydata.w.long <- reshape(mydata.w, direction="long", varying=1:48, v.names="MM")
> 
> ## keep only obs with non-missing values for MM
> mydata.b.long <- na.omit(mydata.b.long)
> mydata.w.long <- na.omit(mydata.w.long)
> 
> 
> lm.b.out <- lm(MM ~ 1, data=mydata.b.long)
> lm.w.out <- lm(MM ~ 1, data=mydata.w.long)
> 
> tab.b <- coef_test(lm.b.out, vcov="CR2", cluster=mydata.b.long$id,
+                    test="Satterthwaite")
> 
> tab.w <- coef_test(lm.w.out, vcov="CR2", cluster=mydata.w.long$id,
+                    test="Satterthwaite")
> WminusB.diff <- tab.w[1,1] - tab.b[1,1]
> WminusB.var <- tab.w[1,2]^2 + tab.b[1,2]^2
> WminusB.SE <- sqrt(WminusB.var)
> WminusB.z <- WminusB.diff / WminusB.SE
> WminusB.pval <- 2*(1-pnorm(abs(WminusB.z)))
> 
> cat("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
> cat("                Averaging Over Only Black Photos\n")
> cat("----------------------------------------------------------------\n")
> cat("Part 1:  E[y | i %in% I^W] - E[y | i %in% I^B] > 0\n")
> cat("\n")
> cat("----------------------------------------------------------------\n")
> cat("Avg. MM score among white respondents\n")
> print(as.data.frame(tab.w), digits=4)
> cat("\nAvg. MM score among black respondents\n")
> print(as.data.frame(tab.b), digits=4)
> cat("\n\n----------------------------------------------------------------\n")
> cat("WminusB.diff     WminusB.SE     WminusB.z    p-value (2-tailed)\n")      
> cat("----------------------------------------------------------------\n")
> cat(WminusB.diff, "       ", WminusB.SE, "    ", WminusB.z, "       ",
+     WminusB.pval, "\n")
> cat("----------------------------------------------------------------\n")
> cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n\n")
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> ## Part 2
> ## E[y | i %in% I^W_{M1}] - E[y | i %in% I^B_{M1}] >
> ##    E[y | i %in% I^W_{M2}] > E[y | i %in% I^B_{M2}]
> ## where the expectation is over all respondents assigned to M1 or M2,
> ## all positions in S, and all photos in J^B
> 
> mydata <- read.csv("../Data/Cleaned\ Data/ScaleRaceSpring2017clean.csv")
> 
> 
> ## just black respondents
> mydata.b <- mydata[mydata$R.race == "Black",]
> mydata.b <- mydata.b[!is.na(mydata.b$R.race),]
> ## just white respondents
> mydata.w <- mydata[mydata$R.race == "White",]
> mydata.w <- mydata.w[!is.na(mydata.w$R.race),]
> 
> 
> ## keep all photos but nothing else
> inds <- 1:24
> keep.vars.template <- c(paste("X", inds, "_Q329", sep=""),
+                         paste("X", inds, "a.id", sep=""))
> keep.vars.notemplate <- c(paste("X", inds, "_Q339", sep=""),
+                           paste("X", inds, "a.id", sep=""))
> 
> mydata.b.m1 <- mydata.b[, keep.vars.template]
> mydata.b.m2 <- mydata.b[, keep.vars.notemplate]
> mydata.w.m1 <- mydata.w[, keep.vars.template]
> mydata.w.m2 <- mydata.w[, keep.vars.notemplate]
> 
> ## convert factor to character
> for (i in 25:48){
+     mydata.b.m1[,i] <- as.character(mydata.b.m1[,i])
+     mydata.w.m1[,i] <- as.character(mydata.w.m1[,i])    
+     mydata.b.m2[,i] <- as.character(mydata.b.m2[,i])
+     mydata.w.m2[,i] <- as.character(mydata.w.m2[,i])    
+ }
> ## get rid of white photos
> for (i in 1:24){
+     i2 <- i+24
+     indic.b.m1 <- grepl("mw", mydata.b.m1[,i2])
+     indic.w.m1 <- grepl("mw", mydata.w.m1[,i2])
+     indic.b.m2 <- grepl("mw", mydata.b.m2[,i2])
+     indic.w.m2 <- grepl("mw", mydata.w.m2[,i2])
+     mydata.b.m1[indic.b.m1, i] <- NA
+     mydata.w.m1[indic.w.m1, i] <- NA    
+     mydata.b.m2[indic.b.m2, i] <- NA
+     mydata.w.m2[indic.w.m2, i] <- NA    
+ }
> mydata.b.m1 <- mydata.b.m1[,1:24]
> mydata.w.m1 <- mydata.w.m1[,1:24]
> mydata.b.m2 <- mydata.b.m2[,1:24]
> mydata.w.m2 <- mydata.w.m2[,1:24]
> 
> 
> 
> ## reshape into long format
> mydata.b.m1.long <- reshape(mydata.b.m1, direction="long",
+                             varying=1:24, v.names="MM")
> mydata.b.m2.long <- reshape(mydata.b.m2, direction="long",
+                             varying=1:24, v.names="MM")
> mydata.w.m1.long <- reshape(mydata.w.m1, direction="long",
+                             varying=1:24, v.names="MM")
> mydata.w.m2.long <- reshape(mydata.w.m2, direction="long",
+                             varying=1:24, v.names="MM")
> 
> ## keep only obs with non-missing values for MM
> mydata.b.m1.long <- na.omit(mydata.b.m1.long)
> mydata.b.m2.long <- na.omit(mydata.b.m2.long)
> mydata.w.m1.long <- na.omit(mydata.w.m1.long)
> mydata.w.m2.long <- na.omit(mydata.w.m2.long)
> 
> lm.b.m1.out <- lm(MM ~ 1, data=mydata.b.m1.long)
> lm.b.m2.out <- lm(MM ~ 1, data=mydata.b.m2.long)
> lm.w.m1.out <- lm(MM ~ 1, data=mydata.w.m1.long)
> lm.w.m2.out <- lm(MM ~ 1, data=mydata.w.m2.long)
> 
> tab.b.m1 <- coef_test(lm.b.m1.out, vcov="CR2", cluster=mydata.b.m1.long$id,
+                    test="Satterthwaite")
> tab.b.m2 <- coef_test(lm.b.m2.out, vcov="CR2", cluster=mydata.b.m2.long$id,
+                    test="Satterthwaite")
> 
> tab.w.m1 <- coef_test(lm.w.m1.out, vcov="CR2", cluster=mydata.w.m1.long$id,
+                    test="Satterthwaite")
> tab.w.m2 <- coef_test(lm.w.m2.out, vcov="CR2", cluster=mydata.w.m2.long$id,
+                    test="Satterthwaite")
> 
> WminusB.m1.diff <- tab.w.m1[1,1] - tab.b.m1[1,1]
> WminusB.m1.var <- tab.w.m1[1,2]^2 + tab.b.m1[1,2]^2
> WminusB.m1.SE <- sqrt(WminusB.m1.var)
> WminusB.m1.z <- WminusB.m1.diff / WminusB.m1.SE
> WminusB.m1.pval <- 2*(1-pnorm(abs(WminusB.m1.z)))
> 
> WminusB.m2.diff <- tab.w.m2[1,1] - tab.b.m2[1,1]
> WminusB.m2.var <- tab.w.m2[1,2]^2 + tab.b.m2[1,2]^2
> WminusB.m2.SE <- sqrt(WminusB.m2.var)
> WminusB.m2.z <- WminusB.m2.diff / WminusB.m2.SE
> WminusB.m2.pval <- 2*(1-pnorm(abs(WminusB.m2.z)))
> 
> WminusB.diff.diff <- WminusB.m1.diff - WminusB.m2.diff
> WminusB.diff.diff.var <- WminusB.m1.var + WminusB.m2.var
> WminusB.diff.diff.SE <- sqrt(WminusB.diff.diff.var)
> WminusB.diff.diff.z <- WminusB.diff.diff / WminusB.diff.diff.SE
> WminusB.diff.diff.pval <- 2*(1-pnorm(abs(WminusB.diff.diff.z)))
> 
> 
> cat("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
> cat("                Averaging Over Only Black Photos\n")
> cat("----------------------------------------------------------------\n")
> 
> cat("Part 2: (E [y | i in I^W_M1 ]  -  E [y | i in I^B_M1 ]) -\n")
> cat("         (E[y | i in I^W_M2 ]  -  E [y | i in I^B_M2 ]) > 0\n") 
> cat("\n")
> cat("----------------------------------------------------------------\n")
> cat("Avg. MM score among white respondents & no template\n")
> print(as.data.frame(tab.w.m1), digits=4)
> cat("\nAvg. MM score among black respondents & no template\n")
> print(as.data.frame(tab.b.m1), digits=4)
> cat("\n\n")
> cat("Avg. MM score among white respondents & template\n")
> print(as.data.frame(tab.w.m2), digits=4)
> cat("\nAvg. MM score among black respondents & template\n")
> print(as.data.frame(tab.b.m2), digits=4)
> 
> cat("\n\n----------------------------------------------------------------\n")
> cat("WminusB.m1.diff  WminusB.m1.SE  WminusB.m1.z   p-value (2-tailed)\n")      
> cat("----------------------------------------------------------------\n")
> cat(WminusB.m1.diff, "       ", WminusB.m1.SE, "    ", WminusB.m1.z, "       ",
+     WminusB.m1.pval, "\n")
> cat("----------------------------------------------------------------\n")
> 
> cat("\n\n----------------------------------------------------------------\n")
> cat("WminusB.m2.diff  WminusB.m2.SE  WminusB.m2.z   p-value (2-tailed)\n")      
> cat("----------------------------------------------------------------\n")
> cat(WminusB.m2.diff, "       ", WminusB.m2.SE, "    ", WminusB.m2.z, "       ",
+     WminusB.m2.pval, "\n")
> cat("----------------------------------------------------------------\n")
> 
> cat("\n\n****************************************************************\n")
> 
> cat("----------------------------------------------------------------\n")
> cat("WminusB.dif.dif  WminusB.dif.dif  WminusB.dif.dif.z   p-value (2-tailed)\n")      
> cat("----------------------------------------------------------------\n")
> cat(WminusB.diff.diff, "       ", WminusB.diff.diff.SE, "    ", WminusB.diff.diff.z, "       ",
+     WminusB.diff.diff.pval, "\n")
> cat("----------------------------------------------------------------\n")
> 
> cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n\n")
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> sink()
> 
> 
> proc.time()
   user  system elapsed 
  1.484   0.105   1.580 
